* FigureA1.do  09/25/22
* Display cross-sections of Rx share and opioid death rates

* get aggregate cross sections
use "PopSeriesAnnual.dta" if race=="All" & gender=="All" & age=="0+" & year==2020, clear
replace year = 2021
tempfile fi21
save `fi21', replace
use "PopSeriesAnnual.dta" if race=="All" & gender=="All" & age=="0+" & year<=2020, clear
append using `fi21'
sort gender age race state place year
save `fi21', replace

use "DeathSeriesAnnual.dta" if race=="All" & gender=="All" & age=="0+" & ucd=="Drug" & (mcd=="T402T403" | mcd=="T400T401T404" | mcd=="Opioid"), clear
sort gender age race state place year
merge gender age race state place year using `fi21'
drop if race!="All"
tab _merge
drop _merge

* aggregate to subperiods
gen byte y20112021 = year>2010
collapse (mean) deaths pop, by(mcd state y20112021)

* reshape mcd codes
encode mcdtitle, gen(mcdnum)
reshape groups mcdnum 1 2 3
reshape cons state y20112021 pop
reshape var deaths
reshape wide
ren deaths1 deathsall  // will less than Rx + Im due to poly drug use
ren deaths2 deathsim
ren deaths3 deathsrx
gen polysh = (deathsim + deathsrx)/deathsall - 1
gen deathrtall = deathsall*100000/pop
gen imsh = deathsim/(deathsim+deathsrx)
gen rxsh = 1 - imsh
gen deathrtrx = (1-imsh)*deathrtall
gen deathrtim = imsh*deathrtall
gen lndeathrtall = log(deathrtall)

gen str18 stshort="New England" if state=="Division 1: New England"
replace stshort="Middle Atlantic" if state=="Division 2: Middle Atlantic"
replace stshort="East North Central" if state=="Division 3: East North Central"
replace stshort="West North Central" if state=="Division 4: West North Central"
replace stshort="South Atlantic" if state=="Division 5: South Atlantic"
replace stshort="East South Central" if state=="Division 6: East South Central"
replace stshort="West South Central" if state=="Division 7: West South Central"
replace stshort="Mountain" if state=="Division 8: Mountain"
replace stshort="Pacific" if state=="Division 9: Pacific"

#delimit ;
scatter rxsh deathrtall if y20112021==0, mlabel(stshort) mlabsize(vsmall) xscale(log) xlabel(3 5 10 20 30) || scatter rxsh deathrtall if y20112021==1, mlabel(stshort) mlabsize(vsmall) xscale(log)
  title("Figure A1.  The cross-area relation between opioid death" "rates and their Rx-Im composition, by subperiod") subtitle(" ")
  legend(order(1 "1999-2010" 2 "2011-2021") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("Opioid fatality rate per 100K, log scale") ytitle("Rx share of opioid deaths")
  note("Notes: Each subperiod is an average of the corresponding years." "Each year's data is the aggregate of all demographic groups in the Census division." "Source: CDC Wonder.", size(vsmall));
#delimit cr

